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ABSTRACT 

We detect and quantify significant numerical biases in the determination of 
the slope of power laws with Salpeter (or similar) indices from uniformly-binned 
data using chi-square minimization. The biases are caused by the correlation 
between the number of stars per bin and the assigned weights and are especially 
important when the number of stars per bin is small. This result implies the 
existence of systematic errors in the values of IMFs calculated in this way. We 
propose as an alternative using variable-size bins and dividing the stars evenly 
among them. Such variable-size bins yield very small biases that are only weakly 
dependent on the number of stars per bin. Furthermore, we show that they 
allow for the calculation of reliable IMFs with only a small total number of stars. 
Therefore, they are a preferred alternative to the standard uniform-size binning. 

Subject headings: methods: numerical — methods: statistical — stars: luminos- 
ity function, mass function 

1. Introduction 

The data explosion of the last decade created by the growth of telescope diameters, 
number of space missions, and detector surface area has produced a comparable growth in 
the amount of objects that can be observed with a single exposure or analyzed in a single 
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paper. In the field of stellar astronomy, this has translated into thousands or larger amounts 
of stars being included in a single color-magnitude diagram or a single luminosity or mass 
function. Indeed, many articles aim at measuring the initial mass function (or IMF) of a 
cluster or of a stellar population by measuring the luminosity of each star, converting the 
results into masses (using isochrone fitting or directly counting the number of stars between 
evolutionary tracks), and fitting a power law or similar function to the masses organized in 
bins using chi-square minimization. 

A known numerical effect takes place when data are binned and a function is fitted to 
the outcome: a bias in the derived parameters can be (and usually is) present if there is a 
low number of objects in some of the bins (see, e.g. Bevington & Robinson 1992; Nousek 
& Shue 1989). The problem originates in the strong anticorrelation between the data and 
the weights in chi-square minimization for data with Poisson uncertainties (Wheaton et al. 
1995, the same should be true for binomial uncertainties, which in most cases can be well 
approximated by Poisson ones, see e.g. Bevington & Robinson 1992). The bias can be 
especially large when some bins have or 1 counts; in the former case some algorithms 
set the weights to correspond to that of 1 count, but that does not eliminate the effect. 
Furthermore, the bias can also be present when the number of objects per bin is not too 
small if some bins are more heavily populated than others (Wheaton et al. 1995). The 
problems associated with this bias are usually correctly dealt with in some astronomical 
fields, such as in high-energy astronomy, but it is surprising to find that many of the articles 
on the IMF appear to ignore them. Thus, it is not strange to find articles in which an IMF 
is calculated with some bins having just a few stars and other bins with tens or hundreds 
of them (see Humphreys & McElroy 1984; Massey et al. 1989; Bejar et al. 2001 for some 
examples). Such a binning scheme introduces biases in the calculated IMF slopes (Kroupa 
2001; Elmegreen 2004). 

There are different approaches to minimizing or eliminating binning biases. Kearns 
et al. (1995) suggest using the fitted number of counts (instead of the real number of counts) 
to calculate the weight of each bin, a method that is iterative by nature since one does not 
know the fitted number of counts a priori. An alternative recommended by D'Agostino & 
Stephens (1986) involves two measures: 

• Define bins of variable size in x and adjust them in such a way that each one of them 
has approximately the same number of objects. The reasoning behind this idea is to 
assign the same statistical weight to each bin and thus minimize biases. 

• If A is the total number of objects, then divide the data into m 2 • A 2 / 5 bins. 

In this article we set out to quantify the importance of binning biases for the deter- 
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mination of the IMF by means of a series of simple numerical experiments. We point out 
that our results should be relevant not only to IMF calculations but to similar problems as 
well, such as the calculation of cluster mass or luminosity distributions. We will start by 
using "standard" uniform-size bins and then we will explore the use of non-uniform ones, as 
suggested by D'Agostino & Stephens (1986). 



2. Definitions and experiments 

Following the Scalo (1986) notation, we define the IMF for stars f(m) such that f(m)-dm 
is the number of stars formed at the same time in some volume of space with mass in the 
interval m to m + dm. Throughout this article, we will assume that the true IMF follows a 
power-law distribution: 



dn 
dm 



f{m) = A ■ m 7 , 



(1) 



where the Salpeter slope is given by 7 = —2.35 (Salpeter 1955). 



The number of stars created in the interval [m a , m b ] can be obtained by integrating 
Eqn. 1: 



r m b dri f mb 
/ - — dm = A ■ m? 'dm 
Jm n dm J„ 



(2) 



to obtain: 



A 



n\Z h a = n(m b ) - n(m a ) = —j—^ ■ [m 6 7+1 - m a 7+1 ] ; 7 -1. 



(3) 



Defining Am, = m b — m a and Xi = ma + mb it is possible to rewrite the last expression as: 



TV- = 77 \ mb 
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(4) 



Finally, we obtain the logarithm of both sides to obtain an expression in the form of yi 
(logarithm of the number of stars in bin i) as a function of Xi (mass at the center of bin i as 
defined in a linear scale): 
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y* = iogio N i = lo gio 



7 + 1 



A 



( 



Xi + 



Ami 



2 




Am 8 



2 



) 




(5) 



For our numerical experiments we used a random number generator to produce 1000 
realizations with 1000 stars each, distributed according to Eqn. 1 with a Salpeter slope 
between m = 6.31 M (log lo m/M = 0.8) and m = 158.49 M (log lo m/M = 2.2). 
The extremes were selected as typical for studies of massive stars (see, e.g. Oey & Clarke 
2005), where the IMF has been measured to be close to the Salpeter value under different 
circumstances, but the conclusions of this work are not expected to change if somewhat 
different values are used. The 1000 realizations were manipulated in two different ways to 
test different conditions: 

• Selecting the first 30, 100, 300, and 1000 (all) stars to simulate different total number 
of stars N = J2 N i- 

• Binning the data into 3, 5, 10, 30, or 50 bins to test whether there is an optimal bin 
size. Given that the D'Agostino & Stephens (1986) recommendation yields 7.8, 12.6, 
19.6, and 31.7 bins for our four values of N, the bins selected here should be able to 
test its validity for our experiments. 

We performed experiments using three types of binning: 

1. Uniform bin size in logarithmic scale (or log 10 [(:Tj + Ami/2)/(xi — Ami/2)] constant) 
with the left edge of the first bin, rad OWI1 , equal to 10 a8 M and the right edge of the 
last bin, m up , equal to 10 22 M (i.e. the input values used for the random generation 
of all realizations). 

2. Approximately constant number of stars in each of the bins with rridown = 10 a8 M 
and m up = 10 2 ' 2 M . This rule cannot be made exact since it is not possible to divide 
e.g. 100 stars in 30 bins with the same number of stars in each. In such a case, we 
design our bins so that they contain e.g. 3, 4, 3, 3, 4, 3. . . stars. 

3. Same as the previous experiment but with m down and m up determined from the data 
in each realization. 

The issue of the weights to be applied to the problem of IMF fitting has generated some 
confusion. Some authors have even decided to skip it altogether by not using any (Massey 



- 5 - 



et al. 1995), but such a strategy can yield large biases. When calculating an IMF we are 
assuming that an underlying physical law determines the true probability distribution f(m), 
which we could measure with arbitrary precision if an infinite number of stars were generated 
from it. In reality, we are limited to analyze finite samples of size N drawn from f(m). It is 
easy to show that, under such circumstances, the value of iVj follows a binomial distribution 
characterized by N and pi, where: 



rXi+Arm/2 ., . , 
_ J Xi - Arm/2 J\ m ) am ( 

Pl P up f{m)dm ' 1 } 

which has a mean of Npi and a variance of Npi(l—pi). Note that the Poisson approximation 
corresponds to the case where pi <C 1; however, if a bin contains a large fraction of the objects, 
then such approximation is no longer correct. The problem here, as noted by Wheaton et al. 
(1995), is that the true Pi needs to be determined from the true f(m), which is unknown! 
One obvious alternative is to substitute pi by Ni/N, the value derived from the data, so that 
the estimate for the uncertainty associated to iVj becomes y/~Ni(N — N{) /N. That is the 
alternative we will use, keeping in mind that this step is precisely the source of the potential 
biases we have previuosly referred to (Wheaton et al. 1995) and which we are trying to 
minimize by selecting bins with similar numbers of stars. 

With the above considerations and given that (a) = log 10 Ni, and (b) the weight Wi 
associated to a value with uncertainty (i.e. y i ±s i ) 1 when using a chi-square minimization 
algorithm is 1/sf, we have: 

Wl (N-N t )(\og 10 er [) 

which has the consequence of giving zero weight to bins with no stars. This has the advantage 
of yielding no numerical problems there, as opposed to using Ni for the independent variable, 
which produces infinite weights 2 . 

For each of the k (k = 1,1000) realizations in each of the N + bin-size combinations, 
the resulting data (x iy Am;, y h w^) was fitted using a chi-square minimization algorithm 



1 Si is the standard deviation derived from the parent distribution of yi and its measured value, not from 
the expected one. Therefore, what we are describing here is the modified chi-square minimization method 
(Wheaton et al. 1995). 

2 Note, however, that in our case w; becomes infinite when Ni — N but, of course, trying to fit a function 
to a histogram that has all data in a single bin is an ill-defined problem. 
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programmed in IDL in order to obtain 7^ and its uncertainty (the algorithm also yields 
Ak and its uncertainty, which will be ignored here). In order to test for a possible algorithm 
dependence, we calculated our results using both (a) CURVEFIT, which uses gradient expan- 
sion, can be found in the standard IDL distribution, and is originally based on Numerical 
Recipes (Press et al. 1986); and (b) Craig B. Markwardt's 3 MPCURVEFIT, which is based on 
MINPACK-1, available from Netlib. No significant differences were found between the two 
algorithms. 

We analyzed the results for the ensemble of 1000 realizations in each N + bin size 
combination. More specifically, we calculated the mean value for the power-law exponent: 



1000 



k=l 

the mean uncertainty in the power-law exponent: 



1000 



k=l 

and the bias normalized with respect to the uncertainty: 



1000 

& = ^V 7fc + 2 - 35 . (io) 

1000 h °* 

The value of b is the logical criterion to judge the existence of biases. If |6| <C 1, then 
the fitting method will be unbiased because it will yield values that will be larger than the 
real one on ^50% of the occasions and smaller on another m 50%. If, on the other hand, 
|6| ~ 1 or larger, a significant bias will exist. 

Of course, we want to test whether our experiments yield biased or unbiased results in 
order to decide which binning technique is the optimum one. However, a fitting method can 
be unbiased but still yield an incorrect uncertainty estimate. A more complete test would 
be to analyze the distribution of the quantity (7^ + 2.35) /a^, whose mean is given by b and 
whose standard deviation is given by: 



3 See http : //cow . physics . wise . edu/~craigm/idl/idl . html. 
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A binning technique that is both unbiased and yields correct uncertainty estimates 
should produce a distribution for (7^ + 2.35) jo~\. that resembles a Gaussian with 6 = 
and f3 — l. Under such circumstances, we can predict that the true value of the slope of 
the IMF will be e.g. within 7^ — and 7^ approximately 34.1% (or e.g. within 7^ and 
7^ + 2<jfc approximately 47.7%) of the times the experiment is executed. If (3 is significantly 
different from 1 but b 0, then the technique will yield unbiased results with incorrect 
uncertainty estimates. If, on the other hand, (3 ~ 1 but b is significantly different from 
zero, the technique will be biased (e.g. will have a systematic error) but its uncertainty 
estimate (sometimes called the random error) will be correct, as previously mentioned. Of 
course, the worst case scenario implies values of b and (3 significantly different from zero and 
one, respectively, in which case the technique produces both systematic errors and incorrect 
uncertainty estimates. 

Finally, we expect the results from different experiments to depend on the existence of 
bins with few or no stars. In order to analyze that effect we define iVj jm i n to be the mean 
Ni in the bin with the lowest number of counts, which for uniform-size bins and a Salpeter 
power law will be the rightmost one. For our variable-size bin experiments, iVj is the same 
for any bin so there is no need to select an specific one. 



We show in Table 1 the values of 7, a, and b for the first experiment, in which we use 
a uniform bin size. In all of the twenty cases, the value of 7 is found out to be larger that 
—2.350, with the maximum at —1.513. Similarly, b is always positive. In a few cases, the 
bias is small but in most of them it is quite large, b is strongly anticorrelated with N ijmin , as 
we show in Fig. 1. This is an expected behavior, since for bins with a small number of stars 
the difference between the used weight and the real weight (the one derived from the parent, 
not the sample distribution) becomes larger in relative terms. Furthermore, there appears 
to be a critical value around iVj min 10: for lower values biases are quite large while for 
higher ones they are small (though not always negligible). These results confirm what we 
mentioned in the introduction: For a binning scheme where the number of objects per bin is 
highly variable, significant biases can be present even if all bins have more than one star. 



3. Experiment 1: Uniform-size bins 



The existence of significant biases makes the use of uniform-size bins inadvisable for 



Table 1. Results for the uniform bin-size case for 3, 5, 10, 30, and 50 bins. 
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a 










b 
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5 


10 


30 


50 


3 


5 


10 


30 


50 


3 


5 


10 


30 


50 


30 


-2.308 


-2.207 


-2.050 


-1.691 


-1.513 


0.303 


0.275 


0.282 


0.293 


0.296 


0.376 


0.655 


1.181 


2.393 


2.986 


100 


-2.343 


-2.306 


-2.244 


-2.042 


-1.895 


0.155 


0.144 


0.148 


0.154 


0.156 


0.176 


0.376 


0.772 


2.058 


2.988 


300 


-2.345 


-2.335 


-2.316 


-2.231 


-2.157 


0.088 


0.083 


0.085 


0.088 


0.089 


0.121 


0.224 


0.430 


1.384 


2.200 


1000 


-2.345 


-2.344 


-2.339 


-2.314 


-2.289 


0.048 


0.045 


0.046 


0.048 


0.049 


0.151 


0.163 


0.260 


0.766 


1.275 



oo 

I 



- 9- 



3- 5 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 



3.0 



2.5 



2.0 — 



■o 1.5 — 



1.0 



0.5 



0.0 — 



* * 



★ * 



* * 



_Q_ 5 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 



log N 

3 1 i,min 



Fig. 1. — Bias as a function of iVj im i n for the first experiment. Note that Ni jUl i n can be 
smaller than 1 because it is a property derived from the parent distribution. 
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calculating IMFs. However, if no alternative is available, one should prefer those results 
based on a small number of bins, since in that case biases are smaller. Increasing the number 
of bins reduces the uncertainty estimate slightly but at the cost of introducing significantly 
larger biases (i.e. it produces smaller random uncertainties but larger systematic errors). 
This is another consequence of the anticorrelation between iVj imin and b. 

We also performed a check on the values of (3 and found them to be close to 1 in all cases. 
Furthermore, a look at the sample case in Fig. 2 shows that the distribution of (7^ + 2.35) /a k 
is well characterized by a Gaussian with a mean equal to b and a standard deviation of 1.0. 
As indicated in the previous section, this means that bins with uniform size yield correct 
uncertainty estimates (i.e. the problem lies in the systematic, not in the random errors). 

We should point out that the existence of a bias towards a flattening of the IMF for 
small samples when using uniform-size bins was previously detected by Kroupa (2001), who 
called it a sampling bias. As we will see in the next sections, it is possible to get rid of such 
a bias almost completely. 

4. Experiment 2: variable-size bins with evenly divided number of stars per 

bin 

For our second experiment we adopt the first recommendation of D'Agostino & Stephens 
(1986) and use a variable size for our bins designed in such a way as to have a similar number 
of stars per bin. We do so for each of our random realizations in the following way (we use 
the case with 100 stars and 5 bins as an example): [a] fixing m down and m up to be 10 a8 M Q 
and 10 2 - 2 M , respectively; [b] sorting the data so that mi < m 2 < . . . < m 100 ; and [c] fixing 
the limits between bins i and i + 1 to be 0.5 * (m 2 o« + ^2o«+i)- 

We show in Table 2 the values of 7, a, and b for the second experiment. Only nineteen 
cases were used because in one circumstance there were more bins than stars. 7 is very close 
to —2.35 in all cases, with a minimum of —2.403 and a maximum of —2.348. \b\ <C 1 in all 
nineteen cases (maximum value of 0.062), with b being positive in some cases and negative in 
others. Note that the signs of 7 + 2.35 and of b can be different due to the possible existence 
of correlations between the values of 7^ and a k (see also Fig. 5). 

These results indicate that using a variable bin size to include a similar number of stars 
in each bin is a good way of minimizing binning biases. A comparison with the previous 
experiment shows that this is done at no significant cost of increasing a. Regarding the 
second recommendation of D'Agostino & Stephens (1986), we only find a weak dependence 
of b in the number of bins. The robustness of the method is emphasized by the fact that even 
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Fig. 2. — Histogram with the distribution of (7^ + 2.35)/<7fc for the 1000 realizations of the 
first experiment with 300 stars and 50 bins. A Gaussian distribution with mean b = 2.200 
and dispersion of 1.0 is also plotted for comparison. The vertical lines mark the position of 
and of b. 



Table 2. Results for the variable-bin size case for 3, 5, 10, 30, and 50 bins. 
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30 


-2.402 


-2.402 


-2.403 


-2.394 




0.285 


0.285 


0.288 


0.292 




0.012 


-0.024 


-0.046 


-0.029 




100 


-2.365 


-2.367 


-2.366 


-2.364 


-2.366 


0.152 


0.152 


0.154 


0.157 


0.156 


0.015 


-0.018 


-0.024 


-0.029 


-0.030 


300 


-2.353 


-2.352 


-2.353 


-2.355 


-2.354 


0.086 


0.087 


0.088 


0.089 


0.090 


0.034 


0.031 


0.012 


-0.012 


-0.009 


1000 


-2.350 


-2.348 


-2.349 


-2.349 


-2.349 


0.047 


0.047 


0.048 


0.049 


0.049 


0.043 


0.062 


0.045 


0.037 


0.043 
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when 30 stars are divided into 30 bins (i.e. 1 star/bin) no significant biases are detected. 
This result corroborates that the existence of large biases originates in the assignment of 
incorrect weights to each bin and not so much by the fact that one specific bin has a low 
number of stars (Wheaton et al. 1995). 

5. Experiment 3: setting the lower and upper mass limits from the data 

The previous two experiments have an artificial component in them: we are using as 
values for md 0ffI1 and m up the input ones, i.e., the values that are attained only in a sample 
with an infinite number of stars. When using real data, however, those values have to be 
determined. An observer typically fixes the first one using incompleteness criteria (since 
lower-mass stars usually exist but are harder to detect) and the second one is usually un- 
known (and, likely, a quantity one is interested in measuring, see e.g. Oey & Clarke 2005). 
Therefore, we can simulate more realistic conditions by modifying the extremes of our second 
experiment by setting m,j OWI1 — mi — 0.5 * (m 2 — mj, m up = m N + 0.5 * {m N — tojv-i), which 
can be determined directly from the data. 

We show in Table 3 the values of 7, a, and b for the third experiment. Results are 
similar to the ones for the second experiment, with values of 7 between —2.365 and —2.298. 
This is especially so for the cases with 300 and 1000 stars, as expected. Comparing results 
one by one, we find that 7 is always larger here than for the previous experiment, b is now 
always positive and the values of |6| are somewhat higher than in the previous experiment. 
Still, all cases with more than 3 stars per bin have \b\ < 0.134 and even the worst case has 
only |6| = 0.264. In Fig. 3 we see that b is still anticorrelated with Ni, the average number 
of stars per bin. 

These results indicate that fixing the lower and upper mass limits from the data and 
using a variable bin size to include a similar number of stars in each bin is a practical way of 
minimizing binning biases. There is a price in the former of larger biases that has to be paid 
for the lack of knowledge of the ends of the distribution, but that price is small, as one can 
observe by comparing Figs. 1 and 3. We should point out that the second recommendation 
of D'Agostino & Stephens (1986) does not work out for this experiment: if one is interested 
in minimizing biases, then it is preferable to go with a low number of bins (but, as previously 
indicated, biases are never large anyway). 

The values of (3 are found to be close to 1 in all cases and the histogram for the sample 
case in Fig. 4 shows the same behavior as the one in Fig. 2: the distribution of (7fc + 2.35)/(x fc 
is well characterized by a Gaussian with a mean equal to b (and, in this case, close to zero) 



Table 3. Results for the variable-bin size case with data-determined m down and m up for 3, 5, 10, 30, and 50 bins. 
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Fig. 3. — Bias as a function of Ni for the third experiment. Note that the vertical scale for 
the plot is 1/10 that of Fig. 1. 
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Fig. 4. — Histogram with the distribution of (7^ + 2.35) jo^ for the 1000 realizations of the 
third experiment with 300 stars and 50 bins. A Gaussian distribution with mean b = 0.090 
and dispersion of 1.0 is also plotted for comparison. The vertical lines mark the position of 
and of b. 
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and a dispersion of 1.0. Therefore, our technique yields not only a nearly bias-free value for 
the slope of the IMF but also a correct estimate of its uncertainty. 

Another advantage of the method proposed here can be extracted from Table 3: a nearly 
bias-free measurement of the power-law exponent with an uncertainty of less than 0.2 can be 
obtained with only 100 stars. If the number is lowered to 30 stars, then the uncertainty in the 
power-law exponent is close to 0.3. Detailed results for the 30 stars + 5 bins case are shown 
in Fig. 5. Note how the first histogram shows a symmetric distribution while the second 
one is distinctly asymmetric. The difference is explained by the correlation (r = —0.58667) 
between 7& and shown in the bottom plot: chi-square fitting yields larger values of the 
uncertainty in the slope for lower values of the slope itself. Note that this correlation is not 
a problem in itself because the histogram that determines through its mean and dispersion 
whether the technique yields a correct estimate of the IMF slope is the first one, not the 
second. Therefore, we conclude that it is possible to conduct precise studies of the mass 
segregation within a large cluster or to measure the IMF in a small one. 

Different binning schemes can easily yield different values for the same sample. This 
can be seen in Fig. 6, where we show a comparison between the three experiments for a 
single realization. With that in mind, one wonders whether part of the variations in the 
IMF detected by a number of authors (see Elmegreen 2004 for a recent review) are not real 
but simply numerical effects introduced by the different schemes used. In order to test that, 
one would have to reanalyze the data in a uniform manner using an unbiased scheme, such 
as the one presented in this article. 

6. Summary and future work 

We conclude that the binning mechanism proposed in this paper for the fitting of power 
laws with Salpeter slopes yields results that (a) are nearly bias-free and (b) produce correct 
uncertainty estimates, as tested by our numerical simulations. On the other hand, the 
standard uniform-size binning introduces biases that are dependent on the number of stars 
per bin. The power of the technique described here extends to small samples, since we have 
shown that it is possible to obtain accurate values with reasonable precisions for the IMF 
slope even when as few as 30 stars are available for analysis. 

We are finishing an analysis of HST/WFPC2 stellar photometry of the nearby dwarf 
starburst galaxy NGC 4214. In that article we will apply the technique described here in 
order to study the IMF for the massive stars in that galaxy. We will also investigate other 
possible sources of biases in the calculation of the IMF. 




J 1 1 1 1 1 1 1 1 1 1 1 1 II 1 1 1 
_ • .« .. 






K^frV. 



Fig. 5. — Detailed results for the 30 stars + 5 bins case for the third experiment, (top) 
Histogram with the distribution of (7* + 2.35)/<7fc. A Gaussian distribution with mean 
b = 0.134 and dispersion of 1.0 is also plotted for comparison. Vertical lines mark the 
position of and b. (center) Histogram with the distribution of 7^. A vertical line marks 
the position of 7. (bottom) Results for 7^ and at for the 1000 realizations. Lines mark the 
values of 7 and a. 
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Fig. 6. — Comparison between the data and the fitted functions for one of the realizations 
with 30 stars and 10 bins for the three experiments in this article. The top, center, and bot- 
tom panels show the first, second, and third experiments, respectively. Note the differences 
in the size of the error bars for the histograms between the first and the last two experiments. 
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We would also like to point out that, given the purely numerical nature of our analysis, 
our results could be extended to other similar problems. For example, the mass function for 
young stellar clusters can be rather well approximated by a power law with a slope of —2.0 
(see e.g. Fall & Zhang 2001), which is quite close to —2.35, so the same type of biases should 
be present there as well. Note, however, that for a problem where we expect the power 
law to have a radically different exponent (e.g. —10.0), the results in this paper may not 
apply because of the larger disparity between the values of the function at the two extremes. 
In general, we recommend that biases be evaluated for any function fitted to binned data 
through chi-square minimization by means of specific numerical experiments similar to the 
ones in this article. 

We would like to thank the referee, Miguel Cerviho, for very useful suggestions and 
for comments that helped improved this article considerably, especially concerning the issue 
of weights. Support for this work was provided by NASA through grants GO-09419.01-A 
and AR-09553.02-A from the Space Telescope Science Institute, Inc., under NASA contract 
NAS5-26555. 
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